/*-------------------------------------------------------------------------*/
/**
  @file		plotlib_graph.c
  @author	Renaud Wanschoor
  @date	july 2004
  @version	$Revision: 2.3 $
  @brief	contains plib functions that are specific to a graph

  plib functions are used to serialize matchings sequences and results of each
  level of the uncoarsening stage. They are currently used for METIS and PaGrid
  pagridplot also used plib functions.

  there are also utilities used by pagrid in plib, those are error handling,
  memory allocation, priority queues, randum number generator, dll as macros
  io functions, in fact all utils that may used for other programs.

  This files contains functions for ::PlotGraphType management
  
*/
/*--------------------------------------------------------------------------*/

#include "plib.h"

PlotGraphType* PLIBG_Init(void) 
{
  PlotGraphType *graph;

  graph = U_malloc(1,PlotGraphType,"PLIB_InitGraph error code 1\n");
  graph->nvtxs = graph->nedges = -1;
  graph->adjidx = graph->vwgt = graph->adjncy = graph->adjwgt = NULL;
  graph->adjwgtsum = graph->map = graph->match = NULL;
  graph->pro = NULL;
  graph->father = graph->child = NULL;

  MEM_Init(&(graph->mem));
  
  return graph;
}

void PLIBG_Allocate(PlotGraphType* g){
  Mem* mem = &(g->mem);

  MEM_Need(mem, 6*g->nvtxs + 2*g->nedges + 1, int);
  MEM_Allocate(mem);
  
  g->adjidx    = MEM_Request(mem, g->nvtxs + 1, int);
  g->adjncy    = MEM_Request(mem, g->nedges,    int);
  g->vwgt      = MEM_Request(mem, g->nvtxs,     int); 
  g->adjwgt    = MEM_Request(mem, g->nedges,    int); 
  g->adjwgtsum = MEM_Request(mem, g->nvtxs,     int);
  g->map       = MEM_Request(mem, g->nvtxs,     int);
  g->match     = MEM_Request(mem, g->nvtxs,     int);
  g->pro       = MEM_Request(mem, g->nvtxs,     int);

}

void PLIBG_Free(PlotGraphType *graph){
  if(graph!=NULL){
    MEM_Free(&(graph->mem));
    U_free("PLIBG_Free error code 1\n",&graph,LEND);
  }
}

int PLIBG_ReadWhere(PlotGraphType *graph,char* F_output, int id){
  PlotGraphType* ini = PLIBG_Init();
  char out[256];
  FILE* fp;
  int levels = 0;
  if(id == 0)
    sprintf(out,"%s.where",F_output);
  else
    sprintf(out,"%s.where.%d",F_output,id);

  WARNING(((fp = fopen(out,"rb")) == NULL,"MD_PagridWriteWhere error code 2 : cannot open %s\n",out),0);

  while(graph->child != NULL)
    graph = graph->child;

  /* copy last level to have inipartitioning info */
  
  PLIBG_Copy(ini,graph);
  PLIBG_Insert(graph,ini);
  graph = ini;

  do{
    U_unserializeIntArray(graph->nvtxs,fp,graph->pro,LEND);
    graph = graph->father;
    levels++;
  }while(graph != NULL);

  fclose(fp);
  return levels - 1;
}

int PLIBG_Copy(PlotGraphType* dest,PlotGraphType* src){
  int nvtxs = src->nvtxs;
  int nedges = src->nedges;
  //*dest = *src;
  
  dest->wgtflag = src->wgtflag;
  dest->ncon    = src->ncon;
  dest->nvtxs   = src->nvtxs;
  dest->nedges  = src->nedges;
 
  PLIBG_Allocate(dest);  
  
  if(src->adjidx != NULL)
    memcpy(dest->adjidx,src->adjidx,(nedges + 1)*sizeof(int));
  if(src->vwgt != NULL)
    memcpy(dest->vwgt,src->vwgt,nvtxs*sizeof(int));
  if(src->adjncy != NULL)
    memcpy(dest->adjncy,src->adjncy,nedges*sizeof(int));
  if(src->adjwgt != NULL)
    memcpy(dest->adjwgt,src->adjwgt,nedges*sizeof(int));
  if(src->adjwgtsum != NULL)
    memcpy(dest->adjwgtsum,src->adjwgtsum,nvtxs*sizeof(int));
  if(src->match != NULL)
    memcpy(dest->match,src->match,nvtxs*sizeof(int));
  if(src->map != NULL)
    memcpy(dest->map,src->map,nvtxs*sizeof(int));
  if(src->pro != NULL)
    memcpy(dest->pro,src->pro,nvtxs*sizeof(int));
  return 1;
}

int PLIBG_SwitchWhere(PlotGraphType *g,int *where,int **old_where){
  //printf("1 %p %p %p %p\n",g->pro,where,old_where,*old_where);fflush(stdout);
  if(old_where != NULL)
    *old_where = g->pro;
  //printf("2 %p %p %p %p\n",g->pro,where,old_where,*old_where);fflush(stdout);
  g->pro = where;
  //printf("3 %p %p %p %p\n",g->pro,where,old_where,*old_where);fflush(stdout);
  return 1;
}

int PLIBG_Insert(PlotGraphType* current,PlotGraphType* new){
  int i;
  
  if(current->child != NULL)
    current->child->father = new;
  new->child = current->child;
  new->father = current;
  current->child = new;
  
  /* build cmap and where */
  for(i=0;i<current->nvtxs;i++){
    current->map[i]=i;
    current->match[i] = i;
  }
  if(new->child != NULL){
    for(i=0;i<current->nvtxs;i++){
      new->pro[i] = new->child->pro[new->map[i]];
    }
  }
  return 1;
}

int PLIBG_DegreeStats(PlotGraphType *g,double *min,double *max,double *avg,double *stddev, double *sum){
  int i;
  int *d = U_mallocset(g->nvtxs,int,0,"MD_PagridDegreeStats error code 1\n");

  for(i=0;i<g->nvtxs;i++)
    d[i] = g->adjidx[i+1] - g->adjidx[i];
  
  U_intStat(d,g->nvtxs,min,NULL,max,NULL,avg,stddev,sum);
  U_free("MD_PagridDegreeStats error code 2\n",&d,LEND);
  
  return 1;
}

int PLIBG_CutEdgeStats(PlotGraphType *g,double *min,double *max,double *avg,double *stddev, double* sum){
  int *c = g->adjwgtsum;
  U_intStat(c,g->nvtxs,min,NULL,max,NULL,avg,stddev,sum);
  *sum/=2;
  return 1;
}

int PLIBG_WeightStats(PlotGraphType *g,double *min,double *max,double *avg,double *stddev, double *sum){
  int *w = g->vwgt;
  U_intStat(w,g->nvtxs,min,NULL,max,NULL,avg,stddev,sum);
  return 1;
}

int PLIBG_EweightStats(PlotGraphType *g,double *min,double *max,double *avg,double *stddev, double *sum){
  int *w = g->adjwgt;
  U_intStat(w,g->nedges,min,NULL,max,NULL,avg,stddev,sum);
  *sum/=2;
  return 1;
}

int* PLIBG_ReadIniWhere(PlotGraphType *g,char* F_output,int id,int *ptr){
  FILE* fp;
  char buf[256];
  int* where;
  while(g->child != NULL)
    g = g->child;

  if(ptr == NULL)
    where = U_malloc(g->nvtxs,int,"MD_PagridReadIniWhere error code 1\n");
  else
    where = ptr;

  sprintf(buf,"%s.where.%d",F_output,id);
  WARNING(((fp = fopen(buf,"rb")) == NULL,"MD_PagridReadIniWhere error code 2 : cannot open %s\n",buf),NULL);
  where = U_unserializeIntArray(g->nvtxs,fp,where,LEND);
  fclose(fp);
  return where;
}

PlotGraphType* PLIBG_Read(FILE *fp){
  PlotGraphType *graph;
  
  int i, k, l, fmt, readew, readvw, ncon, edge, ewgt = 0;
  int *adjidx, *adjncy, *vwgt, *adjwgt, *adjwgtsum;
  char *line, *oldstr, *newstr;
  int sum;		/*auxiliary variable used for calculate adjwgtsum for each vertex*/
  int degree,maxdegree=0;
  
  graph = PLIBG_Init();

  line = U_malloc(MAXLINE+1,char,"MD_PagridReadGraph error code 1\n");

  /* commented lines */
  do { 
    fgets(line, MAXLINE, fp);
  } while (line[0] == '%' && !feof(fp));

  if (feof(fp)) {
    graph->nvtxs = 0;
    U_free("MD_PagridReadGraph error code 1\n",&line,LEND);
    return graph;
  }

  fmt = ncon = 0;  
  sscanf(line, "%d %d %d %d", &(graph->nvtxs), &(graph->nedges), &fmt, &ncon);

  readew = (fmt%10 > 0);	/*determine readew, if the graph has weights associated with the edges, readew = 1*/
  readvw = ((fmt/10)%10 > 0);		/*determine readvw, if the graph has weights associated with the vertices, readvw = 1*/
  if (fmt >= 100) {
    printf("Cannot read this type of file format!");
    exit(0);
  }
  
  graph->wgtflag = 0;
  if (readew)
    graph->wgtflag += 1;
  if (readvw)
    graph->wgtflag += 2;
  
  if (ncon > 0 && !readvw) {
    printf("------------------------------------------------------------------------------\n");
    printf("***  I detected an error in your input file  ***\n\n");
    printf("You specified ncon=%d, but the fmt parameter does not specify vertex weights\n", ncon);
    printf("Make sure that the fmt parameter is set to either 10 or 11.\n");
    printf("------------------------------------------------------------------------------\n");
    exit(0);
  }
  
  graph->nedges *= 2;
  ncon = graph->ncon = (ncon == 0 ? 1 : ncon);
  
  PLIBG_Allocate(graph);
  
  adjidx    = graph->adjidx;
  adjncy    = graph->adjncy;
  vwgt      = graph->vwgt; 
  adjwgt    = graph->adjwgt; 
  adjwgtsum = graph->adjwgtsum;
  
  
  if (! readvw){	/*If the graph has not weights associated with the vertices, initialize all the weights of vertices to be 1*/
    for (i=0; i<graph->nvtxs; i++)	  graph->vwgt[i] = 1;
  }
  if (! readew){	/*If the graph has not weights associated with the edges, initialize all the weights of edges to be 1*/
    for (i=0; i<graph->nedges; i++)	  graph->adjwgt[i] = 1;
  }
  
  /* Start reading the connectivities of graph file */
  for (adjidx[0]=0, k=0, i=0; i<graph->nvtxs; i++) {
    do {
      fgets(line, MAXLINE, fp);
    } while (line[0] == '%' && !feof(fp));
    sum =0;
    oldstr = line;
    newstr = NULL;

    if (strlen(line) == MAXLINE){ 
      printf("\nBuffer for fgets not big enough!\n");
      exit(0);
    }

    if (readvw) { 	/*If there are more than one weights
			  associated with a vertex, read them*/
      for (l=0; l<ncon; l++) {
	vwgt[i*ncon+l] = (int)strtol(oldstr, &newstr, 10);
	oldstr = newstr;
      }
    }
    
    degree = 0;
    for (;;) {		/*read each edge of vertex i*/
      /*get the label of vertex of a edge*/
      edge = (int)strtol(oldstr, &newstr, 10) - 1;	
      oldstr = newstr;
      
      if (readew) {	/*If there is a weight associated with a 
			  edge, read it*/
	ewgt = (int)strtol(oldstr, &newstr, 10);
	oldstr = newstr;
      }
      
      if (edge < 0)
	break;
      
      adjncy[k] = edge;	/*write the label of vertex of a 
				  edge to adjncy list of edge k*/
      if (readew){	/*If there is a weight associated 
			  with a edge, the adjwgt[k] will 
			  be specified as ewgt, otherwise, 
			  the adjwgt[k] will be assigned as 1*/
	adjwgt[k] = ewgt;
	sum += ewgt;	/*adding adjwgtsum*/
      }
      else{
	sum += 1;	/*adding adjwgtsum*/
      }
      k++;		/*increase the label of vertex with 1*/
      degree++;
    } 
    if(degree>maxdegree)
      maxdegree = degree;
    adjwgtsum[i] = sum;	/*assign value to adjwgtsum[i]*/
    adjidx[i+1] = k;		/*specify the end of index of adjacent vertices. Recall that the adjacent vertices of vertex i is storing in adjidx[i] to adjidx[i+1]*/
  }
  
  if (k != graph->nedges) {
    printf("------------------------------------------------------------------------------\n");
    printf("***  I detected an error in your input file  ***\n\n");
    printf("In the first line of the file, you specified that the graph contained\n%d edges. However, I only found %d edges in the file.\n", graph->nedges/2, k/2);
    if (2*k == graph->nedges) {
      printf("\n *> I detected that you specified twice the number of edges that you have in\n");
      printf("    the file. Remember that the number of edges specified in the first line\n");
      printf("    counts each edge between vertices v and u only once.\n\n");
    }
    printf("Please specify the correct number of edges in the first line of the file.\n");
    printf("------------------------------------------------------------------------------\n");
    exit(0);
  }

  U_free("MD_PagridReadGraph error code 2\n",&line,LEND);
  return graph;
}

PlotGraphType *PLIBG_Unserialize(int* tool,
				 char *F_ser,PlotGraphType *graph,
				 char **F_graph,char **F_grid,
				 int* level, 
				 int* CType, int* WType, int *RType,
				 int* CSeed, int *ISeed, int *RSeed){
  int i;
  FILE *fp, *fpgraph;
  PlotGraphType *fgraph;

  fp = PLIB_UnserializeInit(tool, F_ser, F_graph, F_grid,
			    level, CType, WType, RType,
			    CSeed, ISeed, RSeed);

  CALL(*tool,Randini(*CSeed));

  if(graph == NULL){
    WARNING(((fpgraph = fopen(*F_graph,"r")) == NULL,"MD_PagridReadSerial error code 2 : cannot open %s\n",*F_graph),NULL);    
    graph = PLIBG_Read(fpgraph);
    fclose(fpgraph);
  }

  fgraph = graph;

  for(i=0;i<*level;i++){
    U_unserializeIntArray(graph->nvtxs,fp,graph->map,LEND);
    graph = PLIBG_Buildcoarse(*tool,graph);
  }

  U_unserializeIntArray(graph->nvtxs,fp,graph->pro,LEND);

  fclose(fp);

  return graph;
}

PlotGraphType* PLIBG_Buildcoarse(int tool, PlotGraphType* graph){
  int i, idx, j, v, u, k, maxdegree=0;
  
  int fnvtxs     = graph->nvtxs;
  int fnedges    = graph->nedges;
  int* vwgt      = graph->vwgt;
  int* adjidx    = graph->adjidx;
  int* adjncy    = graph->adjncy;
  int* adjwgt    = graph->adjwgt;
  int* adjwgtsum = graph->adjwgtsum;
  int* map       = graph->map;
  int* match     = graph->match;
  int* mapidx    = U_mallocset(fnvtxs,int,-1,"Pagrid_Buildcoarse error code 1\n");

  PlotGraphType *child;
  int *mapadj;
  int *record;

  int 
    nedges,  ncon,
    cnedges, cnvtxs = 0;
  
  int 
    *cadjncy, *cadjwgt, *cadjwgtsum, *cvwgt, *cadjidx;  

  /* rebuild match */
  for(i=0;i<fnvtxs;i++){
    idx = map[i];

    if(idx+1>cnvtxs)
      cnvtxs = idx+1;

    if(mapidx[idx] == -1)
      mapidx[idx] = i;
    else{
      match[i] = mapidx[idx];
      match[mapidx[idx]] = i;
      mapidx[idx] = -2;
    }
  }

  for(i=0;i<fnvtxs;i++){
    idx = map[i];
    if(mapidx[idx] != -2){
      mapidx[idx] = -2;
      match[i] = i;
    }
  }

  /* rebuild guide */
  CALL(tool,Perm(fnvtxs, mapidx));
  
  /* contract */
  child         = PLIBG_Init();
  child->nvtxs  = cnvtxs;
  ncon          = child->ncon  = graph->ncon;
  child->nedges = fnedges;

  PLIBG_Allocate(child);
  
  cadjidx    = child->adjidx;
  cadjncy    = child->adjncy;
  cvwgt      = child->vwgt;
  cadjwgt    = child->adjwgt; 
  cadjwgtsum = child->adjwgtsum;
  
  cnedges = cadjidx[0] = 0;

  /*Initialize the vertices list of the coarser graph (-1 => not contracted)*/
  for (i =0; i<cnvtxs; i++)
    child->vwgt[i] = -1;
  
  record  = U_mallocset(cnvtxs,int,-1,"Contract error code 2\n");
  mapadj  = U_malloc(fnedges,int,"Contract error code 3\n");

  /* map the adjacency matrix of the finer graph */
  for (i = 0; i < fnedges; i++)
    mapadj[i] = map[adjncy[i]];

  /*****************************/
  /* Contraction               */
  /*****************************/
  cnvtxs = 0;
  for (i = 0; i < fnvtxs; i++) {
    v = mapidx[i];
    
    /*if the label to which vertex v is mapped is equal to cnvtxs, 
      collapse it and the vertex to which it matches*/
    if (map[v] != cnvtxs)
      continue;
      
    if (ncon == 1)
      cvwgt[cnvtxs] = vwgt[v];

    /*****************************/
    /* map v in coarser graph    */
    /*****************************/
    cadjwgtsum[cnvtxs] = adjwgtsum[v];	/*Initilize adjwgsum of 
					  vertex cnvtxs in coarser 
					  graph to the adjwgtsum of 
					  vertex v in the father graph*/
    
    nedges = 0;	/*Initialize number of edges of map[v] to be zero*/
    
    /* add all edges of v in the finer graph */
    for (j = adjidx[v]; j < adjidx[v+1]; j++) {
      k = mapadj[j];	
      /* if edges is not already existing add it */
      if (record[k] == -1) {
	cadjncy[nedges] = k;
	cadjwgt[nedges] = adjwgt[j];
	record[k] = nedges++;
      /* else only update its weight */
      } else {
	cadjwgt[record[k]] += adjwgt[j];
      }
    }
    
    /**********************************************/
    /* map u (matched vertex) in coarser graph    */
    /**********************************************/
    u = match[v];	
    if (v != u) {			
      if (ncon == 1)
	cvwgt[cnvtxs] += vwgt[u];	
      
      cadjwgtsum[cnvtxs] += adjwgtsum[u];	
      
      /* add all edges of u in the finer graph */
      for (j = adjidx[u]; j < adjidx[u+1]; j++) {
	k = mapadj[j];
	/* if edges is not already existing add it */
	if (record[k] == -1) {
	  cadjncy[nedges] = k;
	  cadjwgt[nedges] = adjwgt[j];
	  record[k] = nedges++;
	/* else only update its weight */
	} else {
	  cadjwgt[record[k]] += adjwgt[j];
	}
      }
      
      /* Remove the contracted edge if existing */
      if ((j = record[cnvtxs]) != -1) {
	cadjwgtsum[cnvtxs] -= cadjwgt[j];	/* subtract its weight */
	cadjncy[j] = cadjncy[--nedges];		/* replace with last edge */
	cadjwgt[j] = cadjwgt[nedges];		
	record[cnvtxs] = -1;
      }
    }
    /*********************/
    /* clean and update  */
    /*********************/
    for (j=0; j<nedges; j++)
      record[cadjncy[j]] = -1;	
    
    cnedges += nedges;			/*add up the number of edges 
					  found in a iteration*/	
    cadjidx[++cnvtxs] = cnedges;	/*add up the edges found in a 
					  iteration to adjidx of the 
					  coarser graph*/
    
    cadjncy += nedges;		/*move the pointer preparing 
				  for receving new adjacent edges*/
    cadjwgt += nedges;		/*move the pointer preparing 
				  for receving new adjacent edges*/
    if(nedges>maxdegree)
      maxdegree = nedges;
  }
  /************/
  /* finalize */
  /************/
  child->nedges = cnedges;
  
  /*attach child to graph*/
  graph->child = child;
  child->father = graph;
  
  /*free tmp arrays*/
  U_free("Contract error code 3\n",&mapidx,&mapadj,&record,LEND);

  return child;
}

void PLIBG_FreeMultiple(PlotGraphType* graph){
  PlotGraphType* g;
  while(graph->father != NULL)
    graph = graph->father;

  do{
    g = graph;
    graph = graph->child;
    PLIBG_Free(g);
  }while(graph != NULL);
}

double* PLIBG_VtxsDegree(PlotGraphType *g,double *ptr){
  double* d;
  int i;
  if(ptr == NULL)
    d = U_malloc(g->nvtxs+1,double,"MD_PlotGetDegree error code 1\n");
  else
    d = ptr;
  d[0] = g->nvtxs;
  for(i=0;i<g->nvtxs;i++)
    d[i+1] = g->adjidx[i+1] - g->adjidx[i];
  return d;
}

double* PLIBG_VtxsWeight(PlotGraphType *g,double *ptr){
  double* w;
  int i;
  if(ptr == NULL)
    w = U_malloc(g->nvtxs + 1,double,"MD_PlotGetWeight error code 1\n");
  else
    w = ptr;
  w[0] = g->nvtxs;
  for(i=0;i<g->nvtxs;i++)
    w[i+1] = g->vwgt[i];
  return w;
}

double* PLIBG_VtxsCutedge(PlotGraphType *g,double *ptr){
  double* c;
  int i;
  if(ptr == NULL)
    c = U_malloc(g->nvtxs + 1,double,"MD_PlotGetCutEdge error code 1\n");
  else
    c = ptr;
  c[0] = g->nvtxs;
  for(i=0;i<g->nvtxs;i++)
    c[i+1] = g->adjwgtsum[i];
  return c;
}

int* PLIBG_BuildWhere(PlotGraphType *cg,PlotGraphType *fg){
  /* init */
  int 
    i,*tmpwhere,*iniwhere,
    *cwhere = U_malloc(fg->nvtxs,int,"MD_PagridBuildWhere error code 1\n"),
    *where = U_malloc(fg->nvtxs,int,"MD_PagridBuildWhere error code 2\n");
  PlotGraphType *g;

  WARNING((cg->pro == NULL,"MD_PagridBuildWhere error code 3 : no processor information\n"),0);
  iniwhere = fg->pro;
  memcpy(cwhere,cg->pro,cg->nvtxs*sizeof(int));

  /* rock 'n roll */
  do{
    g = cg->father;
    for(i=0;i<g->nvtxs;i++)
      where[i] = cwhere[g->map[i]];
    tmpwhere = cwhere;
    cwhere = where;
    where = tmpwhere;
    cg = g;
  }while(g!=fg);

  /* this is the end */
  g->pro = cwhere;
  U_free("MD_PagridBuildWhere error code 4\n",&where,LEND);  

  return iniwhere;
}

int* PLIBG_PartEdges(PlotGraphType *g,int vtx, int* ptr){
  int* 
    edges,
    i,
    pos = 1;

  if(ptr == NULL)
    edges = U_malloc(g->nedges + 1,int,"MD_PagridPartEdges error code 2\n");
  else
    edges = ptr;

  for(i=g->adjidx[vtx]; i<g->adjidx[vtx+1]; i++){
    if(g->pro[vtx] == g->pro[g->adjncy[i]]){
      edges[pos++] = g->adjncy[i];
    }
  }
  edges[0] = pos - 1;
  return edges;
}

int PLIBG_NbProc(PlotGraphType *g){
  int min = -1;
  int max = 0;
  int i, *where = g->pro, size = g->nvtxs;

  WARNING((g->pro == NULL,"MD_MetisNbProc error code 1 : no processors info\n"),0);
  for(i=0;i<size;i++){
    if(where[i] < min || min == -1)
      min = where[i];
    if(where[i] > max)
      max = where[i];
  }
  return max - min + 1;
}
